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We extend the method of Silvestrelli [P. L. Silvestrelli, J. Chem. Phys. 139, 054106 (2013)] to 
approximate long-range van der Waals interactions at the density functional theory level based on 
maximally localized Wannier functions combined with the quantum harmonic oscillator model, to 
periodic systems. Applying this scheme to study London dispersion forces between graphene and 
water layers, we demonstrate that collective many-body effects beyond simple additive pair-wise 
interactions are essential to accurately describe van der Waals forces. 


I. INTRODUCTION 

Graphene, which consists of a single layer of carbon 
atoms arranged in a sp^ honey-comb structure,has 
been subject of a great number of studies, owing to its 
unique properties that were experimentally observed and 
theoretically predicted. These include, but not limited 
to, electric field^ and quantum Hall effects,'*^ ultra-high 
carrier mobility,® superior thermal conductivity,® high 
mechanical strength,^ electron-hole puddles, and sen¬ 
sitivity to adsorbents.Although, there is a large effort 
focusing on controlling and adapting graphene for future 
use in real devices,some rather fundamental aspects 
still remain incompletely understood. Among others, is 
the interaction between graphene sheets and individual 
water molecules on top of them.^^“^® This interaction can 
cause a doping effect on graphene and change the density 
of states near the Fermi level, which therefore is a rather 
promising effect for potential applications.^^ However, 
weak dispersion interactions play an important role in 
the binding between carbon atoms and water molecules. 
Hence, there seems to be a controversy in the reported 
results, in particular whether or not graphene is trans¬ 
parent with respect to van der Waals (vdW) interactions 
between the substrate beneath graphene and the water 
droplet on top of it.^®“^' Therefore, an accurate treat¬ 
ment of vdW interactions between water and graphene 
is required. The difficulty of conventional local or semi¬ 
local density functional theory (DFT) to quantitatively 
describe these interactions originates from the long-range 
nature of dispersion forces.^® Even though a large num¬ 
ber of different empirical and first-principles techniques 
have been suggested to calculate vdW interactions within 
DFT calculations,there is still great demand for an 
accurate and, at the same time, computationally very 
efficient first-principles technique to account for vdW 
forces in large systems. 

Here, we extend the recently introduced method of 
Silvestrelli®^ which combines Quantum Harmonic Oscil¬ 
lator (QHO) model with Maximally Localized Wannier 


Functions (MLWF),®^ to periodic systems. Using this 
approach, we investigate the vdW interactions between 
an extended layer of water and graphene with respect to 
the number of layers. The question on the additivity of 
vdW interactions between the water and graphene layers 
is discussed in detail. 

This remaining of the paper is organized as follows. 
The theoretical method is outlined in Section If, while 
the computational details are specified in Section HI. The 
eventual results are shown and discussed in Section IV. 


II. METHOD OF CALCULATION 

Considering MLWFs as charge distributions in real- 
space, it is directly possible to compute the correspond¬ 
ing dipole moments.®® Based on this information, the dy¬ 
namical electron-correlation effects that arise from many- 
body instantaneous long-range interactions of the oscil¬ 
lating dipoles can be quantified by the QHO modeU®’^®’®^ 

TV TV N 

i=l i=l i>j = l 

( 1 ) 

where N is the number of MLWFs, while Xi = 
with mi being the masses and Ci the displacements of 
the oscillators from equilibrium. The frequency and po¬ 
larizability of the oscillators are denoted as uJi and at, 
while Tij is the dipole-dipole interaction tensor. Know¬ 
ing the centers and spreads Si of the MLWFs, the re¬ 
spective polarizabilities ai ~ jSf and characteristic fre¬ 
quencies (Mi ~ \/V ai can be calculated, where 7 is a 
proportionality constant and Zi the atomic number.®®’®^ 
Moreover, Tij is modified to allow for orbital overlap at 
short distances, by introducing a short-range damping 
function for the bare Coulomb potential.®® The energies 
of all N 3-dimensional QHOs can be found by diagonal- 




2 


izing a 3N x 3N matrix C that is defined as 

c.. = UJ^ (2a) 

Ci^j — LOiCOj sjO-i Oij Tij , ( 2 ^) 

which contains N'^ 3x3 matrix-blocks corresponding to 
the individual MLWFs. The vdW correction can then be 
obtained by 


Evdw = 


3N 



2=1 


( 3 ) 


where Xi are eigenvalues of the correlated system, while 
uji are the aforementioned characteristic frequencies of 
the dipole moments attributed to the MLWFs. 

Due to the fact that decays relatively quickly with 
respect to the distance between the MLWFs, instead of a 
genuine Ewald sum, the interactions between the MLWFs 
in the unit cell and those in its periodic images are taken 
into account by considering a finite buffer region around 
the unit cell, where the MLWFs of the original unit cell 
are replicated in all spatial directions. To that extend we 
decompose the vdW interaction energy for the extended 
system into a sum of vdW interaction energies of 

the MLWFs in the unit cell , the interaction 

energy between the MLWFs in the unit cell and those in 
the buffer zone E^^^, and interaction energy between 
the MLWFs in the buffer region only, which is denoted 
as E^2w- vdW interaction energy of the extended 
system then reads as 


J^ext _ pUC-UC I pUC-b , j^b-b 
^vdW — ^vdW "I” ^vdW "I” ^vdW^ 


( 4 ) 


from which the last term has to be subtracted to yield the 
desired vdW interaction energy of the original system, i.e. 

T^tot T^ext J7'b—b 

^vdW — ^vdW ^vdW' 

As a consequence, the total energy including the vdW 
correction is E*°* = E*£prp + total 

DFT energy. 


III. COMPUTATIONAL DETAILS 

In the following we have considered two systems, a sin¬ 
gle layer of graphene (SLG), as well as bilayer graphene 
(BLG), both with a 100 molecule water slab on top. The 
graphene, which consisted of 128 carbon atoms per layer, 
was placed in a periodic orthorhombic simulation box 
parallel to the xy-p\a,ne with a large 35 A vacuum por¬ 
tion along the perpendicular z-direction. The MLWF 
centers of the eventual system are shown in Fig. 1 

In order to obtain the atomic configuration of the wa¬ 
ter molecules, the ab-initio molecular dynamics (AIMD) 
simulation were performed in the canonical ensemble 
at 300 K using the second-generation Gar-Parrinello 
method of Kiihne et al.®®’®® as implemented in the Gaus¬ 
sian and plane wave®^ DFT code CP2K/Quickstep.^® 



FIG. 1: The MLWFs centers of the extended system (unit cell 
plus additional buffer region) of the SLG-water system at 3 A 
vertical separation. 



FIG. 2: The vdW energy as a function of buffer region width. 


In AIMD simulation carbon atoms of graphene where 
fixed. The interatomic interactions were described 
by DFT^® employing the Perdew-Burke-Ernzerhof 
(PBE) exchange-correlation functional®^, Goedecker- 
Teter-Hutter pseudopotentials®^’®^ and a double-C Gaus¬ 
sian basis set with one addition al set of polarization 
functions.®® The spread of the Wannier orbitals was min¬ 
imized using the scheme of Berghold et al.®^ 

To estimate the size of the buffer zone in Eig. 2 the 
vdW energy as defined in Eq. 5 is shown as a function 
of additional buffer width, d. We found that d=12 A is 
sufficient to adequately converge the vdW and to capture 
most of the relevant dispersion interactions, while at the 
same time keeping the size of the C matrix manageable. 


IV. RESULTS AND DISCUSSION 


To demonstrate the impact of periodic boundary con¬ 
ditions (PBC), Fig 3(a) displays the vdW energy as a 
function of the distance between the water slab and SLG, 
which is defined as the vertical gap between the near¬ 
est H atom of the water slab and the xy-plane. As can 
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FIG. 3: (a) The vdW interaction energy as a function of dis¬ 
tance between the water slab and the SLG with (blue tri¬ 
angles) and without (red squares) the MLWFs of the buffer 
zone, (b) Relative difference between and 



FIG. 4: Interaction energy between the water slab and 
the SLG as a function of their distance, using PBE with¬ 
out vdW correction (black circles), PBE-D3 (green circles), 
PBE-l-i?^*^^'^ (red squares), and PBE-l-if*^^^^ (blue trian¬ 
gles). 


be seen in Fig. 3(b), the vdW energy with and with¬ 
out our correction for PBC differ by ~34%. In Fig. 4 
the interaction energy Ei^t as obtained from bare PBE, 
PBE-D3, PBE-I-E^^*^'^ and PJiE+El°^yy are shown as 
a function of the vertical distance between the water slab 
and SLG. It is apparent that without any vdW correc¬ 
tion PBE hardly binds, with the associated equilibrium 
distance being about 3.0 A. The corresponding equi¬ 
librium distance for both, PBE-D3 and PBE+El°^y^^ are 
2.5 A and 2.6 A for PBE-|-E^‘^^‘", which is in good 
agreement with previous results obtained by others using 
a large variety of different methods such as polarizable 
and non-polarizable force fields, DFT calculations in¬ 
cluding dispersion corrections, calculations based on the 
random phase approximation, but also with local MP2 


FIG. 5: Interaction energy between the water slab and 
SLG/BLG as a function of distance between the upper 
graphene and the water layer. The PBE results for the SLG 
(BLG) are shown in solid black circles (triangles), while the 
PBE-D3 calculations are denoted by green (orange) spheres 
(squares). The corresponding biding curves in case of the 
present correction for SLG and BLG are depicted 

as red squares and blue triangles, respectively. 


and CCSD(T) calculations of a single water molecule on 
a hydrogen-terminated graphene flake.However, at 
variance to the latter, here a disordered water slab is 
considered, which comprises a large number of different 
water orientations towards the graphene layer. Never¬ 
theless, for most of the water molecules, one 0-H bond 
is preferably pointing towards the hydrophobic surface, 
as has been observed in previous AIMD simulations.®^’®® 
The eventual interaction energy averaged over all water 
molecules is approximately 82 meV per water molecule. 

The results for the interaction energy Eint between 
the water slab and BLG with AB stacking at the exper¬ 
imentally observed separation of 3.34 A,®^^®® are shown 
in Fig. 5 and compared with the values for SLG. Due to 
the fact that considering BLG including the buffer region 
would results in a rather large C matrix that needs to be 
diagonalized, we have confined ourselves to vdW correc¬ 
tion due to MLWFs in the unit cell. As before, employing 
the PBE functional without any vdW correction, water 
and BLG barely binds. More interesting, we find that 
in comparison to the SLG system, Eint is reduced by 
about ~30%. However, in the case of PBE-D3 calcula¬ 
tion, the reduction of Eint due to the second graphene 
layer is just ~5 %, while for E^^^^ the binding ener¬ 
gies between SLG and BLG systems are essentially iden¬ 
tical. Nevertheless, in the latter case, the equilibrium 
distance changes from 2.6 A for SLG to 2.5 A for BLG. 
The PBE-D3 calculations predict an equilibrium distance 
of 2.5 A for both of the considered systems, while for the 
bare PBE functional, the equilibrium distance reduces 
from 3.0 A to 2.9 A for SLG and BLG, respectively. 
In any case, for all the computational methods we have 
considered here, the binding energy between the water 
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slab and SLG is at least at large as for BLG. This is to 
say, that the vdW interactions are non-additive and in 
fact are screened by the additional graphene layer, which 
immediately suggest that the electronic structure of an 
individual sheet is changed dramatically when interact¬ 
ing with other layers. To that extend we have calcu¬ 
lated the z-component of the molecular dipole moment 
of both graphene systems without water using the Berry 
phase operator for periodic systems and the centers of the 


MLWFs.^^’®^ In the case of the SLG, the z-component 
of the dipole is unsurprisingly zero, while in the case of 
BLG, each layer exhibits a dipole moment of ^3.65 D, 
though in opposite directions. We conclude by noting 
that the latter is a manifestation that the whole is more 
than the sum of its constituents, which highlights the im¬ 
portance to explicitly consider the electronic structure of 
the full interacting system to embrace the subtle many- 
body effects of the vdW interaction. 
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